Numerical evolution of multiple black holes with accurate initial data 
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We present numerical evolutions of three equal-mass black holes using the moving puncture ap- 
proach. We calculate puncture initial data for three black holes solving the constraint equations 
by means of a high-order multigrid elliptic solver. Using these initial data, we show the results for 
three black hole evolutions with sixth-order waveform convergence. We compare results obtained 
with the BAM and AMSS-NCKU codes with previous results. The approximate analytic solution 
to the Hamiltonian constraint used in previous simulations of three black holes leads to different 
dynamics and waveforms. We present some numerical experiments showing the evolution of four 
black holes and the resulting gravitational waveform. 
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I. INTRODUCTION 

The gravitational n-body problem is an old and im- 
portant problem which dates back to 1687, when Isaac 
Newton's "Principia" was published. The special case for 
n = 3 was studied by Euler, Lagrange, Laplace, Poincare, 
among others (see, e.g., However, solutions of the 

three-body problem have shown a rich complexity and 
are far from being completely understood. From the 
point of view of celestial mechanics, the three-body prob- 
lem is related to the irnportant question of the stability 
of the solar system 0tJ|- In globular clusters n-body 
interactions appear to be important in the formation of 
intermediate-mass black holes [5l-[TTI|. Finally, covering 
large scales in cosmology numerical solutions of the n- 
body problem are used to simulate formation of structure 

In order to solve the three-body problem, scientists 
have developed mathematical tools, and with the de- 
velopment of computers since the 1950's, also numerical 
techniques. On the one hand, there are analytical solu- 
tions for special cases of the three-body problem, for ex- 
ample those due to Euler and Lagrange (see e.g. [ll.[T6|). 
on the other hand, there are solutions which exhibit a 
chaotic behavior. At the beginning of the 20th century 
Sundman found a convergent series solution to the three- 
body problem However, the rate of convergence 
of the series which he had derived is extremely slow, and 
it is not useful for practical purposes. From the point of 
view of dynamical systems, the three-body problem was 
a key system which allowed Poincare to identify many of 
the novel ideas related to the dynamical system theory 
and chaos fl8l |. 

Using post-Newtonian techniques (PN), it is now pos- 
sible to describe the dynamics of n compact objects, up 
to 3.5 PN order (see e.g. For binary systems the 

ADM Hamiltonian has been specialized up to 3.5 PN or- 
der (2^ . and for three bodies there are explicit formulas 
up to 2 PN order j23l - [25l |. Periodic solutions, also known 
as choreographic solutions, were studied using these tech- 



niques [25M27f ■ as well as estimates of the gravitational 
radiation for binary-single interactions [1, H^, [2^ . 

The first complete simulations using general- 
relativistic numerical evolutions of three black holes 
were presented in [s^, HH (see [32| - [3^ for very limited 
early examples of multiple black hole simulations). The 
recent simulations show that the dynamics of three 
compact objects display a qualitative different behavior 
than the Newtonian dynamics. Gravitational waves 
are an extra component in the three-body problem of 
compact objects which enrich the phenomenology of 
the system. The changes in the energy and momentum 
resulting from the gravitational radiation produce a 
difference in the dynamics of the system. There are open 
questions related to the general-relativistic dynamics 
of n compact objects, for example the possible chaotic 
behavior of the dynamics of n black holes, the inverse 
problem in gravitational wave emission, the existence of 
quasi-stationary solutions and their stability, etc. 

In some regard simulations of three or more black holes 
are more sensitive to small changes in the data than bi- 
nary simulations, as also noted in e.g. [sol. IsTj. In a typi- 
cal binary, changing the initial momentum slightly leads 
to a correspondingly small change in the eccentricity of 
the binary. In a black hole triple, changing the initial mo- 
mentum may change the dynamics completely since the 
first black hole can merge with the second or the third 
depending on the initial momentum. Hence, moving up 
from two to three bodies introduces a new feature, but 
of course from the point of view of n-body simulations 
and chaotic systems this is no surprise. 

In this paper we examine the sensitivity of the fully 
relativistic evolutions of three black holes to changes in 
the initial data. We present simulations of three and four 
black holes. The examples for three black holes are some 
of the simpler cases already considered in [sO, |3l|. A 
more detailed analysis about a possible chaotic behavior 
of the three body problem in general relativity is beyond 
the scope of this work, but would certainly be of interest. 

In (sol , [sij , initial data is specificed using an analytic 
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approximation, which introduces a finite error that does 
not converge to zero with numerical resolution. The rea- 
son to use such initial data is that, although accurate 
initial data for two black holes is readily available, this 
is not the case for more than two black holes. Below we 
show that solving the constraints numerically to obtain 
initial data for an arbitrary number of black holes, the 
result of the evolutions can change dramatically. The ac- 
tual difference between the analytic approximation and 
the numerical initial data is not large (depending on the 
initial parameters), but, as expected, even small differ- 
ences can lead to large changes for multiple black hole 
orbits. 

The paper is organized as follows. In Sec. we re- 
view the puncture method [SS^-'i^, which is the basic 
approach that we use to solve the initial data problem. 
This is followed by a description of our new code, Ol- 
LIPTIC, designed to solve the constraint equations of the 
3-1-1 formalism numerically. Olliptic implements a par- 
allel multigrid algorithm on nested regular grids, with 
up to eighth order finite differencing. In Sec. Ill CI we 
present our results for three test cases, for the initial 
data of a single puncture, and for two and three punc- 
tures. The evolution of three black holes is presented 
in Sec. IIIIl where we compare BAM results directly with 
the AMSS-NCKU code ^ and indirectly with the pre- 
vious results of [30l Isij . Finally, we perform simulations 
of four black holes to show that the same techniques work 
for more than three back holes (our largest simulation in- 
volves 42 black holes 37]). We conclude with a discussion 
in Sec. llVl 



II. INITIAL DATA 

Under a 3-1-1 decomposition, the Einstein equations 
split into a set of evolution equations and constraint 
equations, namely the Hamiltonian and momentum con- 
straints (see, e.g., [42| - |43 | for reviews). In vacuum the 
constraint equations read as follows: 



(1) 
(2) 



Following the conformal transverse-traceless decompo- 
sition approach, we make the following assumptions for 
the metric and the extrinsic curvature: 



K,. 



It3 = ^tln 



(3) 
(4) 



where A^^ is trace free. We choose an initially flat back- 
ground metric, = J^j, and a maximal slice, if = 0. 
The last choice decouples the constraint equations (HJ-© 
which take the form 



0, 



(5) 
(6) 



Bowen and York [46| have obtained a non-trivial so- 
lution of Eq. ([5|) in a Cartesian coordinate system (2;*), 
which by linearity of the momentum constraint can be 
superposed for any number of black holes (here the in- 
dex n is a label for each puncture): 
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where r„ := + yn + is the Levi-Civita tensor 

associated with the flat metric, and Pi and Si are the 
ADM linear and angular momentum, respectively. 

The Hamiltonian constraint ^ becomes an elliptic 
equation for the conformal factor. The solution is spli t 
as a sum of a singular term and a finite correction u [40] , 



(8) 



with M — !■ as r„ — > 00. The function u is determined by 
an elliptic equation on M.^ and is C°° everywhere except 
at the punctures, where it is C^. The parameter to„ is 
called the bare mass of the 71th puncture. 



Numerical Method 



where R is the Ricci scalar, Kij is the extrinsic curvature 
and K its trace, 7^^- is the 3-metric, and Vj the covariant 
derivative associated with 7ij . 



A. Puncture method 

The constraints can be solved, for example, with the 
puncture method of ji^l. N black holes are modeled by 
adopting the Brill-Lindquist wormhole topology [45| with 
N+1 asymptotically flat ends which are compactified and 
identified with points on R^. The coordinate singular- 
ities at the points resulting from compactification are 
referred to as punctures. 



In order to solve Eq. (|6]) numerically, we have written 
Olliptic, a parallel computational code to solve three 
dimensional systems of non-linear elliptic equations with 
a 2nd, 4th, 6th, and 8th order finite difference multigrid 
method. The elliptic solver uses vertex-centered sten- 
cils and box-based mesh refinement that we describe be- 
low. We use a standard multigrid method |47 H51J with a 
Gauss-Seidel Newton relaxation algorithm (e.g. |52l|V 

The numerical domain is represented by a hierarchy of 
nested Cartesian grids. The hierarchy consists of L -|- G 
levels of refinement indexed byZ = 0,...,L-|-G — 1. A 
refinement level consists of one or more Cartesian grids 
with constant grid-spacing hi on level I. A refinement 
factor of two is used such that hi = /iG'/2l'~'^l . The grids 
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are properly nested in that the coordinate extent of any 
grid at level I > G is completely covered by the grids at 
level l — l. The level / = G is the "external box" where the 
physical boundary is defined. We use grids with I < G 
to implement the multigrid method beyond level I = G. 

The parallelization approach that we use is block de- 
composition, in which each domain is divided into rect- 
angular regions among the processors such that the com- 
putational work load is balanced. For levels I > G every 
domain uses p/2 buffer points at the boundary of the do- 
main (here p indicates the order of the finite difference 
stencil). Levels with I < G contain a single point at the 
boundary. For every face of the three dimensional rect- 
angular domain we use these points for different purposes 
(see Fig.[l]): 

1. If the face is on the outside of the global domain, we 
use the points as a refinement boundary (or phys- 
ical boundary if Z = G); the boundary conditions 
are explained below. 

2. If the face is in the internal part of the global do- 
main, then we use ghost zones of the neighbor- 
ing processors to update information of the buffer 
points. 

3. If the face is defined with symmetry, we use a reflec- 
tion condition to calculate the values at the bound- 
ary. 

Olliptic can be used with three symmetries: oc- 
tant {—x,—y,—z) — >■ (a;,?/, z), quadrant {—x,—y,z) 
{x,y,z) and bitant {x,y,—z) {x,y,z). We use the 
negative part of the domain to define the computational 
grid, because that increases the performance of the re- 
laxation method somewhat since the resulting order of 
point traversal helps propagating boundary information 
into the grid. 

For the "physical" or outer boundary we require that 
u ^ A as r ^ oo. The standard condition used in this 
case is an inverse power fall-off. 



u{r) = A + 



B 



for r > 1, g > 0, 



(9) 



where the factor B is unknown. It is possible to get an 
equivalent condition which does not contain B by calcu- 
lating the derivative of (O with respect to r, solving the 
equation for B and making a substitution in the original 
equation. The result is a Robin boundary condition: 



r duix) 

u{x) + ^ 

q or 



A. 



(10) 



The implementation of the boundary condition was a 
key point to get accurate solutions, so we describe our im- 
plementation in some detail. Rather than taking deriva- 
tives in the radial direction as is required by ([TU|) . we 
take derivatives only in the direction normal to the faces 
of our rectangular domain. At the edges of the bound- 
ary, we use a linear combination of the derivatives along 
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FIG. 1: 2D representation of a domain divided between 2 
processors and with bitant symmetry. Face (A) is internal to 
the global domain and handles a reflection symmetry. Face 
(B) is internal to the domain and manages the communication 
with the second processor. Faces (C) and (D) are physical 
boundaries, where we impose a Robin boundary condition in 
the normal direction n (see text). 



the normals of the two adjacent faces. At the corners, 
we use a linear combination of the derivatives for the 
three adjacent faces. In the computation, we first apply 
the boundary condition to the interior of the boundary 
faces, then compute derivatives inside the faces to up- 
date the edges, and then compute derivatives inside the 
edges to obtain boundary data at the corners. We use a 
one sided finite difference stencil of order p and a Newton 
iteration method to update the values on the boundary. 
For example, for the face (D) (see Fig. [Ij, the equations 
are 



dRijk ■— 



^■IJk 



^ijk 



ijk 



n+1 n Rijk 



ijk 



(11) 

(12) 
(13) 



where -D^^ is the forward difference operator of order p 
in the y-direction, Uijk, rijk and yijk are the values of u, 
r, and y, respectively, at the lattice location (i, j, fe), and 
n is an iteration index. For example, in the case p = 2 
we obtain 



4Mi,j + l.t+Mi,j + 2,fc 

2Ay 



-3/2Ay. 



(14) 
(15) 



Note that Eq. (fTTj) is linear in M^fc, so in fact the al- 
gorithm is equivalent to implementing the explicit finite 
difference formula, with the advantage that its imple- 
mentation is easier. Since there are p/2 boundary buffer 
points, we have to specify a method to obtain more than 
one buffer point. In our implementation the method is 
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TABLE I: Grid setups used for tests and puncture initial data. 
Ir := I — G is the number of inner refinement levels, L is the 
length of the numerical domain and hmin is the grid size in 
the finest level (see text for details about each system). 



System 


Levels 


Length 


Grid size 






L 


hmin 


Test 1 


1-4 


4.8 * 2'"^ 


1/20 


Test 2 & 3 


3 


20.0 


{0.1,0.09, ...,0.02} 


1-puncture 


5 


40.0 


{5/256, 1/64,5/384} 


2-punctures 


7 


40.0 


{1/16,1/32,1/64} 


3-punctures 


7 


50.0 


5/64 



TABLE IL Results of test 1, where p is the order of the stencil 
which we use to solve the equation, I is the number of refine- 
ment levels, R is the relative error calculated in the finest 
level and AR is the comparison with the reference solution. 



p 


2 


4 


6 




8 




I 


R AR 


R AR 


R 


AR 


R 


AR 






(xlO-^) 


(xlO 


-10) 


(xlO 




1 


3.83 


4.36 


9.29 




3.05 




2 


4.57 0.74 


12.56 8.20 


32.81 


23.52 


67.25 


64.20 


3 


5.23 1.40 


15.42 11.06 


105.37 


96.08 


207.75 


204.70 


4 


5.54 1.71 


16.93 12.56 


139.53 


130.24 


284.75 


281.70 



stable if we update the values of the boundary points 
from the inside to the outside of the domain. First, in- 
side points are used to get the first boundary point using 
the one-sided derivative. Then the stencil is shifted by 
one from the inside to the outside, including the first 
boundary point to compute data at the second boundary 
point, and so forth. 



C. Results 

1. Analytic test problems 

We test Olliptic with three simple elliptic equations 
using the following procedure. Given the solution U'^ on 
a mesh with grid-spacing h and an elliptic operator C^, 
we calculate a source p'^ which satisfies the equation 



(16) 



and then we solve the equation to obtain numerically. 
In this way it is possible to calculate the error 



(17) 



where | • | is a suitable norm. We summarize the grid 
setup for our tests and puncture initial data in Table HI 
The goals of the first test were to estimate the error in- 
troduced by the refinement method and to investigate the 
effectiveness of the algorithm to solve non-linear equa- 
tions. We have solved the equation 



V'^U{x) + U{xf = pi{x) 
U{x) = ee-5^-^ 



for 
for 



xen, (18) 

xedn, (19) 



where is the three-dimensional Laplace operator, and 
n is the interior of a rectangular domain. The solution 
given is a Gaussian function with amplitude e — 0.004, 
in this case we use a Dirichlet boundary condition. We 
have solved the equation with a single level of refinement 
in a cube of length L =4.8, and with mesh size dx = 
dy = dz = 0.05. Using this solution as reference, we 
solve Eq. (|18p . increasing the number of levels up to 3 



external boxes. Due to the Dirichlet boundary condition 
the numerical solution is exact at the boundary. We use 
the norm Loo to calculate the relative error. 



R 



(20) 



and as measurement of the error introduced by the re- 
finement method, we calculate the difference between the 
error using more than one refinement level and the refer- 
ence solution, AR = \R{1 > 1) - i?(/ = 0)|. The results 
are summarized in Table |IT1 

The results for the non-linear Eq. show that using 
high order schemes gives a significant improvement in the 
accuracy of the solution. Increasing the order from p to 
p + 2 decreases R by almost three orders of magnitude. 

In order to test the implementation of the Robin 
boundary condition, we use a second trial function. 



V^Uia:)=p2{x) 



U{x) 



tanh(r) 



for f G 17, (21) 
for X e dVL, (22) 



where r :— \x\. The solution [/ is a function which has 
the asymptotic behaviour given by Eq. ^ with A = 0, 
B = 1, and q = I- In this case we look at the convergence 
of our numerical data using 3 levels of refinement in a 
cubic domain of length 20, and using 9 resolutions going 
from 0.1 to 0.02 in the finest level. For a finite difference 
implementation of order p, for h <^ I, we expect 



ChP, 



(23) 



where E'^, is given by Eq. (jl7p using the L2 norm, h 
is the mesh size, and C is constant with respect to h. 
After calculating the logarithm of Eq. (^5)) we get a linear 
function of p, 



ln{E^) ^p\n{h) + C'. 



(24) 



Using this expression with our data and doing a linear 
regression analysis, we estimate the convergence order V 
for our numerical experiment (in the best case V — > p 
as h — 5- 0). As measurement of the error we use the 
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TABLE III: Convergence test for the Robin boundary condi- 
tion. Here p is the order of the finite difference, V, a, and 
are the mean, the standard deviation, and the coefficient of 
variation of the convergence order for our numerical experi- 
ments, respectively, and AV is the relative deviation of our 
results with respect to p. 



P 



V 



AV 



2 2.002 0.0002 0.009% 0.10% 

4 3.994 0.0005 0.013% 0.15% 

6 5.985 0.0013 0.022% 0.26% 

8 7.969 0.0020 0.026% 0.39% 



TABLE IV: Convergence for a solution which is C2°. Here k 
is the exponent given in 



p V a Cy AV 


V (J cy AV 


2 2.003 0.0003 0.013% 0.16% 
4 3.782 0.0088 0.233% 5.46% 
6 3.848 0.0067 0.175% 35.86% 
8 3.836 0.0038 0.098% 52.05% 


1.999 0.0001 0.005% 0.06% 
3.995 0.0003 0.007% 0.12% 
5.715 0.0124 0.216% 4.74% 
5.868 0.0294 0.500% 26.65% 



standard deviation and the coefRcient of variation of our 
data. The resuhs are displayed in Table IIIII 

We have obtained an accurate implementation of the 
boundaries for problem (I?n) - (j22l) . where the difference 
between the theoretical convergence order and the ex- 
perimental one is less than 0.5%. However, note that the 
convergence at the boundary depends on specific proper- 
ties of the test problem. 

For the last analytic test, we verify the accuracy of the 
method for a function which is C2° := C°°(M3 \ {0}). 
The problem to solve was 



V'^U{x) ^ p:i{x) for f 6 ^2, (25) 
U{x) = for X £ dn, (26) 

where we set fc = 3 or fc = 5, r := \x\, and U is 
everywhere except at the origin, where it is C''^^. 

We use the procedure of the second test to estimate the 
convergence order, changing the equation and the bound- 
aries (in this case we use a Dirichlet boundary condition) . 
The result of our numerical experiments (detailed in Ta- 
ble |IV| shows that the overall convergence of the numeri- 
cal solution calculated using a standard finite differencing 
scheme is restricted by the differentiability of the analyt- 
ical solution. The convergence order close to the origin 
(within a few grid points) is the same as the order of 
differentiability and improves significantly moving away 
from the origin. The accuracy of the solution behaves in 
the same manner. 
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FIG. 2: Regular part u of the conformal factor along the Y- 
axis of a single puncture with vanishing spin parameter and 
with linear momentum Py = 0.2M. Shown is a convergence 
test without scaling (left) and with scaling (right) for second- 
order convergence using cf2 = 1.8409. 



2. Single puncture initial data 

After calibrating our code, we calculate the Hamilto- 
nian constraint for a single puncture. We tested the con- 
vergence of our second-order implementation for a single 
boosted puncture (P' = 0.2 (5pf, S' = 0) by looking 
at the value of the regular part u of the conformal fac- 
tor along the F-axis for a cubic domain of length 40M, 
5 levels of refinement, and 3 resolutions hi = (5/8)M, 
/i2 = Ahi/5, and — 2hi/3 in the coarse level. In 
Fig. [21 we show rescaled and unsealed data for positive 
and negative values of Y, respectively. 

We plot the values of \u'^^ - u'^^\ and - for 
y < on the left, and on the right values for F > with 
_u'»3| rnultiplied by a factor cf2 = 1.8409 which cor- 
responds to the proper scaling of second order. The lines 
in the right panel of the plot coincide almost everywhere, 
indicating second order convergence. We also show de- 
tails of a region close to the puncture in the insets. 

We perform a similar test calculating spinning black 
hole initial data (P* = 0, = 0.2(5?,M). Fig. H shows 
the result of the convergence test for this case where we 
found second order convergence again. 

As an example of a high order solution, in Fig. 3] we 
show the convergence test for the eighth order scheme of 
the boosted puncture. In this case the plot shows a drop 
of the convergence ratio close to the puncture. However, 
far from the puncture the convergence behavior is better. 

In Fig. [5] we plot the results for the spinning black 
hole, obtained by using our fourth order implementation. 
Compared to the boosted puncture, in this case we see 
better behavior close to the puncture (the solution of 
the 8th order spinning puncture is similar to the boosted 
case). Far from the puncture the convergence ratio is 
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FIG. 3; Regular part u of the conformal factor along the Y- 
axis of a single puncture with vanishing linear momentum and 
with spin Sy = 0.2M. Shown is a convergence test without 
scaling (left) and with scaling (right) for second-order conver- 
gence using cf2 = 1.8409. 



FIG. 5: Regular part u of the conformal factor along the Y- 
axis of a single puncture with vanishing linear momentum and 
with spin Sy — 0.2AI. Convergence test without scaling (left) 
and with scaling (right) for fourth-order convergence using 
cf4 = 2.7840. 
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FIG. 4: Regular part u of the conformal factor along the Y- 
axis of a single puncture with vanishing spin parameter and 
with linear momentum Py — 0.2M. Eighth-order convergence 
of u is obtained far from the puncture (cfg = 6.4637). 



approximately second order. 

As we saw in our third test and in our numerical 
experiment for a single boosted or spinning puncture, 
the convergence rate of the high order finite differencing 
scheme for functions C2° drops near to 0. This is a well 
known property of high order finite difference schemes 
(e.g. [4^, [531). We review some basics of this effect in 
Appendix |21 Nevertheless, as we show in IIIIBI and in 
the two-punctures test (see below), the numerical solu- 
tion produced by our high order implementation seems 
to be accurate enough to perform numerical evolutions 
of multiple black holes. The errors close to the puncture 
do not modify significantly the convergence during the 
evolution. 



3. Two-puncture initial data 

As a test for a binary system we set the parameters for 
two punctures to xi = -X2 = 3M, ^ -P^ = 0.2 8\M. 
This configuration was studied before using a single- 
domain spectral method [13] . We compared the result of 
our new code with the solution produced by the spectral 
solver. For the spectral solution we use ua — ub — 40 
and — 20 collocation-points (see reference for details 
about the definition of spectral coordinates {A,B,(p)). 
We calculate the multigrid solution in a cubic domain of 
length 40M, 7 levels of refinement and 3 resolutions of 
hi = (1/16)M, /i2 = hi/2 and = /ii/4 in the finest 
level. 

Fig.[6]is a plot similar to Fig. 5 of [5^ . We compare the 
spectral solution with the eighth order multigrid solution. 
The fact that the four lines coincide on the scale of the 
plot (3 resolutions of multigird and one spectral solution) 
indicates that the two methods agree with each other on 
the whole domain. 

Using the same setting we solve the Hamiltonian con- 
straint with the second, fourth, and sixth order stencil 
of the multigrid code. Then we use the highly accurate 
solution of the spectral code as reference to compare with 
the different orders. As we showed before in the case of 
a single puncture, the accuracy close to the puncture de- 
creases. However, the comparison with the spectral code 
(see Fig. [T]) shows that using high order finite differencing 
stencils improves the accuracy of the solution. 

Spectral methods produce in general more accurate so- 
lutions to elliptic equations than those obtained by finite 
difference methods |55[. However, in order to take full 



advantage of the spectral method for punctures, it is nec- 
essary to construct a special set of coordinates. Indeed, 
there exist coordinates in which the conformal correction 
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FIG. 6: Comparison between the numerical solution of the 
Hamiltonian constraint calculated using a single-domain spec- 
tral method and the high-order multigrid solver. The plot 
shows u along the X-axis produced by the spectral code (de- 
noted by us) and using three resolutions calculated with the 
eighth order implementation of the multigrid code (labels 



FIG. 8; Plot of u along X-axis for system 3BH102, compar- 
ing the approximate solution with a second order numerical 
solution. The higher order numerical solutions would not be 
distinguishable from second order in this plot, compare Fig.[7l 



4- Three-puncture initial data 




X(M) 

FIG. 7: Absolute value of the differences between the numer- 
ical solution of u for the second, fourth, sixth, and eighth 
order finite difference implementation and the spectral solu- 
tion. The upper part shows the result for second and fourth 
order, and the bottom part for fourth, sixth, and eighth order. 



u is smooth at the puncture [5J|. Although these coor- 
dinates are in principal applicable for both spectral and 
finite difl:erencing methods, the resulting grids are spe- 
cific to two black holes. Generalizing that approach to 
more than two punctures is an interesting but non-trivial 
challenge that we do not pursue in this work. 

Using finite difference multigrid methods with Carte- 
sian coordinates, one advantage of the puncture construc- 
tion is that it is possible to produce accurate solutions of 
the Hamiltonian constraint for multiple black holes with 
minimal changes to a code prepared for binaries. 



In previous work on the numerical evolution of three 
black holes [sO, Hfj , the Hamiltonian constraint has been 
specified using an approximate solution (see [31 , 56- 59] ) . 
We compare our numerical solution with the approximate 
solution (which we implemented as well) for the set of 
parameters labeled 3BH102 given in Table I of [s^l, see 
our Table E 

In Fig. [5] we show a plot of the solution obtained using 
a cubic domain of length 50M, a mesh size h = 0.5M 
in the coarse level and 9 levels of refinement. The ap- 
proximate solution was calculated in the same numerical 
grid. The result shows a significant difference between 
the two methods, and, as we will show later in llllBl that 
fact leads to a quantitative and qualitative difference for 
evolutions. 



III. NUMERICAL EVOLUTION OF THREE 
BLACK HOLES 

In the mid 1960's, Hahn and Lindquist started the nu- 
merical investigation of colliding black holes [lOl- After 
more than f orty ye ars and a series of breakthroughs start- 
ing in 2005 [6lJ-|65| , the numerical relativity community is 
now able to produce stable black hole inspiral simulations 
and to compute gravitational waves signals. The most 
common formulation used to perform numerical evolu- 
tions of black holes is based on the work of Shibata and 
Nakamura j66i |. and Baumgarte and Shapiro G?] and is 
known as the BSSN formulation. 
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A. Techniques 

We have performed the three black hole simulations 
using the BAM code as described in @, HI], and with 
the AMSS-NCKU code In BAM we use a sixth 
order discretization for the spatial derivatives [69| and 
fourth order accurate integration in time. Initial data are 
provided by the Olliptic code. Gravitational waves are 
calculated in the form of the Newman-Penrose scalar ^4 
according to the procedure described in Sec. Ill of ^6§\. 
We use the BSSN system together with the 1 + log and 
gamma freezing coordinate gauges [70l - l72j | as described 
in "gs"] (choosing in particular the parameter r] = 2/M 
in the gamma freezing shift condition). All the runs are 
carried out with the symmetry (x,y,z) {x,y,—z) in 
order to reduce the computational cost. The Courant 
factor, C := At/hi, seems to be an important ingredient 
to reach convergence. For long evolutions (evolution time 
t > 200), we set C = 1/4, in other cases we use C = 1/2. 

The AMSS-NCKU code is an extended version of 
the code described in [41]. Instead of GrACE, we con- 
structed our own driver combining GH — h and FORTRAN 
90 to implement moving box style mesh refinement. Re- 
garding the numerical scheme dealing with the interface 
of neighbor levels, we closely follow the methods de- 
scribed in [HItI]. AMSS-NGKU can implement both 
the 6 point buffer zone method (68| and interpolation at 
each sub-Rung-Kutta step [t^ ■ Our tests show little dif- 
ference between these two methods. For simplicity, all 
simulations presented here use the 6 point buffer zone 
method. In order to do 3rd order interpolation in time, 
we need three time levels of data. At the beginning of 
the numerical evolution, we use a 4th order Rung-Kutta 
method to evolve the initial data backward one step to 
get the data on time level t — to — At for every mesh 
level [zl]. Here At is different for different mesh levels. 
We have reproduced the results published in [4l| with 
this driver. All of the tests involved fixed mesh refine- 
ment and agree very well with the results obtained with 
GrAGE. As to the Einstein equation solver, we replaced 
the ICN method used in [4l[ by a 4th order Runge-Kutta 
method. The Sommerfeld boundary condition is imple- 
mented with 5th order interpolation. 



TABLE V: Initial data parameters 



Parameter 


3BH1 


3BH102 


FBHSR 


FBH3w 


xi/M 


-2.4085600 


-3.5223800 


-10.0000 


-10.000000 


yi/M 


2.2341300 


2.5850900 


0.0000 


0.0000000 


pf/M 


-0.0460284 


0.0782693 


-0.0100 


0.0000000 


p\/M 


-0.0126181 


-0.0433529 


0.0750 


0.0524275 


nii/M 


0.3152690 


0.3175780 


0.2500 


0.2500000 


X2/M 


-2.4085600 


-3.5246200 


-6.0000 


-6.0000000 


2/2 /M 


-2.1053400 


-2.5850900 


0.0000 


0.0000000 


pIIM 


0.1307260 


-0.0782693 


0.0000 


0.0000000 


pI/m 


-0.0126181 


-0.0433529 


-0.0750 


-0.0524275 


7712 /M 


0.3152690 


0.3175780 


0.2500 


0.2500000 




4.8735000 


7.0447600 


6.0000 


8.0000000 


ya/M 


0.0643941 


0.0000000 


0.0000 


0.0000000 


Pi/M 


-0.0846974 


0.0000000 


0.0000 


0.0000000 


pI/m 


0.0252361 


0.0867057 


0.0750 


-0.0781250 


TTlz 1 iVl 


U.oloZDyU 


U.oloOOOU 


U.zoUU 


u.zouuuuu 


X4,/M 






10.0000 


10.0000000 


yi/M 






0.0000 


0.0000000 


Pl/M 






0.0100 


0.0000000 


pI/m 






-0.0750 


0.0781250 


rm/M 






0.2500 


0.2500000 




-6e-07 - 



60 80 
Time (M) 



B. Results 

With the BAM code we simulate three black holes with 
initial parameters as given in Table [V] In the first exper- 
iments, we focus on runs that use the initial data param- 
eters of runs "SBHl" and "3BH102" in [s^. We evolve 
this data with both the numerical initial data and the 
approximate solution to the conformal factor. We com- 
pare the puncture tracks and the extracted wave forms 
with those produced by the AMSS-NGKU code. The 
puncture tracks give a convenient measure of the black 
hole motion. It is much more cumbersome to compute 
the event horizon, which we do for a simple black hole 



FIG. 9: Real part of ^'4 (mode Z = m = 2) calculated at 
r = 40M for system 3BH1. The lower panel shows the con- 
vergence test for 6th order (cfg = 1.9542). 



triple in [75|- Finally, we discuss some results for the 
evolution of four black holes which show some additional 
properties of multiple black holes evolutions. 



1. Three black holes 

System 3BH1 is a short simulation which is useful 
for convergence tests. We use our sixth order imple- 
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FIG. 10: Puncture tracks for system 3BH102 using approxi- 
mate initial data, (see text, Sec. Ill C"4l) . 



mentation to calculate initial data, a cubic domain of 
length 1052M, 10 levels of refinement and three resolu- 
tion hi = (125/12)M, /i2 = 6hi/7 and = 2hi/3 on 
the finer level. We have obtained roughly sixth order 
convergence for the gravitational waveform, as shown in 
Fig. ini Our results show a \['4 waveform similar to that 
shown in Fig. 16 of [3l'|. For this evolution, we did not 
find a significant difference when using approximate ini- 
tial data or solving the constraint equations numerically. 
However, the accuracy of the numerical initial data is 
important for the numerical result to converge. Our first 
test using BAM's elliptic solver (which is a second-order 
multigrid solver) showed that for long-time simulation it 
is important to improve the accuracy of the initial data 
for multiple black hole evolutions. 

Our second example is black hole configuration 
3BH102, which we consider first for approximate initial 
data, and later for the numerical solution. This set of 
parameters is a system which, starting with approximate 
initial data, leads to trajectories forming a nice figure 
similar to the Greek letters 7, a and t (see Fig. [TUl com- 
puted with BAM). Our convergence test for this sys- 
tem shows sixth-order (see Fig. [T5|) . with small deviations 
from second and fourth order which are consistent with 
the accuracy of the evolution method of our code. 

Comparing with Fig. 3 of [131 ^ there is a small but not- 
icable difference in the puncture tracks of roughly up to 
IM in the coordinates compared to our results. There are 
several possible explanations for this difference. Evolu- 
tions of multiple black holes are sensitive to small changes 
in the grid setup and initial data. We tested possible 
sources of errors, for example introduced by numerical 
dissipation or finite resolution. Changing these lead to 
negligible changes in the trajectories on the scale of the 
plot and do not seem to explain the existing difference. 
However, since the deviation from (s^l does not change 
the qualitative shape of the tracks, we conclude that we 
have consistently reproduced that simulation. 




-8 -4 4 8 

X(M) 

FIG. 11: Puncture tracks for system 3BH102 using ap- 
proximate initial data comparing results for the BAM and 
AMSS-NCKU codes. The difference in the trajectories is 
small, and the results agree in the general shape. Note that 
AMSS-NCKU uses fourth-order spatial discretization instead 
of sixth-order which is implemented in BAM. 
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merical solution of the Hamiltonian constraint with the 8th 
order multigrid method. There is a drastic change in the 
puncture tracks compared to the evolution of the approxi- 
mate initial data, in particular the black holes merge in a 
different order (compare Fig. llOjl . 



Alternatively, we can compare the paths of the punc- 
tures obtained with BAM with those produced by the 
AMSS-NCKU code. The implementation of the ap- 
proximate initial data was done independently for the 
two codes, and in both cases the formula from [sil] is 
used. We see in Fig. [TT] that the results from the two 
codes agree within a maximum difference of about 0.2M 
in the given coordinates, or 2% with respect to an orbital 
scale of lOM. An analysis of the I = m = 2 mode of ^'4 
showed that there are differences in the phase of about 
0.4% and of about 2% in the amplitude. 
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l=2,m=2,r=50M 




50 100 150 200 250 300 350 400 
Time (M) 



FIG. 13: Real part of (mode I — m — 2) calculated at 
r = 50AI for system 3BH102 using approximate initial data. 
The upper panel shows the waveform for 3 resolutions, 
the bottom plot shows sixth-order convergence scaling with a 
factor cfe = 1.9542. 



l=2,ra=2,r=50M 
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Time (M) 



FIG. 14: Real part of (mode I — m — 2) calculated at 
r = 50M for system 3BH102 using the numerical solution of 
the Hamiltonian constraint. The upper panel shows the 
waveform for 3 resolutions, the bottom plot shows sixth-order 
convergence scaling with a factor cfe = 1.9542. Near the first 
merger the order of convergence is closer to fourth-order. 



When comparing codes, recall that the BAM evolu- 
tions use 6th order spatial differencing, AMSS-NCKU 
4th order, and also 4th order for the figures, pointing 
out that there was little difference to an 8th order run. 
Our conclusion is that differences due to resolution are 
small, and they are significantly smaller than the changes 
introduced by replacing the approximate initial data by 
a numerical solution of the Hamiltonian constraint. 

We now focus on the evolution of system 3BH102 
solving the Hamiltonian constraint with the eighth-order 
multigrid method. The ^"4 waveform and convergence 
are shown in Fig. [TH Note that we see approximately 



sixth-order convergence in the waveform except for the 
first merger where the convergence is close to 4th order. 

As shown in Sec. Ill C 4l for system 3BH102 the numeri- 
cal solution of the Hamiltonian constraint differs from the 
approximate prescription. As a consequence, the trajec- 
tories and waveform change. We show the paths followed 
by the punctures for this case in Fig. [121 Instead of the 
grazing collisions of the previous evolution, in this case 
the black holes with labels 2 and 3 merge after a small 
inspiral, producing a higher amplitude in the wave. The 
second merger is almost a head-on collision, which gen- 
erates a smaller amplitude in the wave. Notice that the 
order in which the black holes merge differs from the pre- 
vious evolution. 

Looking at the wave forms, for the approximate initial 
data Fig. [T^] shows a relatively large burst of "junk"- 
radiation which does not converge. Solving the Hamil- 
tonian constraint we see a better convergence behav- 
ior, see Fig. [M] Moreover, the difference in the junk- 
radiation between resolutions using the approximate ini- 
tial data is one order of magnitude bigger than solving 
the Hamiltonian constraint numerically (compare the in- 
sets in Figs. fT3l and fT4)) . 

In the case of a binary system it is possible to produce 
the same evolution for numerical and approximate initial 
data by adjusting the mass parameter [3l|. In the case 
of three black holes, there does not seem to be a simple 
procedure to fit the initial parameters in order to repro- 
duce the same trajectory with both types of initial data. 
We tried changes in the momentum, the mass, and the 
momentum and mass together, looking at the maximum 
of the regular part u of the conformal factor in order to 
reduce the difference between the analytical prescription 
and the numerical data. The result is not satisfactory, 
i.e. we did not find a way to change the parameters of 
the approximate data to better approximate the solution 
of the Hamiltonian constraint, and the large differences 
in the puncture tracks could not be removed. 



2. Four black holes 

Evolution of more than three black holes is possible 
using the same approach. We performed several tests for 
evolutions of multiple black holes. Here we present two 
particular cases using four black holes. 

The first case is the system that we call FBHSR (see 
Table |V] for details about the initial parameters) . We 
start with four equal-mass black holes aligned on the X- 
axis, and with an arbitrary selection of the initial momen- 
tum symmetric with respect to the YZ-plane. The black 
holes follow a rotationally symmetric path (see Fig. [T5]) . 
From the waveform we can distinguish two mergers. The 
first merger is between the black holes labeled BHl and 
BH4, the black hole generated by that merger stays in the 
origin until a triple merger occurs with the black holes 
BH2 and BH3. The second merger is almost a triple 
head-on collision with a small amplitude wave form. 
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FIG. 15: Symmetric configuration of four blaclc lioles. Top; 
Patlis followed by the punctures with a rotational symmetry 
of 180° with respect to the Z-axis. Bottom: Real part of the 
I = m = 2 mode of r»I'4 extracted at r = 50M, the waveform 
shows two mergers (around t = 150M and t = 200M), first 
between the black holes labeled BHI and BH4, and the second 
merger between the resulting black hole and the black holes 
BH2 and BH3. 



i 




100 150 
Time(M) 



250 



FIG. 16: Triple merger of four black holes. Top: Paths fol- 
lowed by the punctures. Bottom: Real part of the I = m = 2 
and / = 2, m = modes of r*&4 extracted at r = 50M. The 
waveform shows three mergers (around t — 55M, t — 95M, 
and t = 170M), which are more easily identified by looking 
at the mode I = 2, m = 0. 



The second case (FBH3w in Table |V]) consists of a 
quick merger of two binary systems (see Fig. I16p . The 
remaining black holes merge in an almost head-on colli- 
sion. Looking at the waveforms (I = m = 2 and Z = 2, 
771 = modes of Re(r^'4)) we can identify the three col- 
lisions, the second with higher amplitude. The initial 
parameters were chosen by trial and error in order to 
produce this kind of waveform. With a larger separation 
and more careful choice of the initial parameters it would 



be possible to find a stronger and more distinctive merger 
signal. However, the above examples are only intended 
to be an illustration of the kind of waveforms that can 
be generated by mergers of four black holes. 



IV. DISCUSSION 

We have presented a numerically elliptic solver, Ol- 
LIPTIC. As a first application, we solve the Hamiltonian 
constraint to obtain numerical initial data for multiple 
black hole evolutions. Olliptic implements a high-order 
multigrid method, which is parallelized and uses a box- 
based mesh refinement. The tests and first applications 
of the code showed that the new code seems to be suf- 
ficiently accurate for our purposes. However, we found 
that close to the puncture the convergence rate is less 
than that desired, which is expected for puncture data 
(see Appendix E]) ■ The drop in the convergence close to 
the punctures is not reflected in the convergence of the 
evolution. Nevertheless, we are considering to modify the 
numerical scheme in order to improve the accuracy close 
to the puncture. 

Wc have shown evolutions of three and four black holes 
which use as initial data solutions to the Hamiltonian 
constraint generated with the new elliptic solver. We 
compare with results for a certain analytic approximation 
for the initial data. In the case of three black holes, 
the dynamics resulting from approximate data is different 
from the dynamics produced by evolutions which satisfy 
the Hamiltonian constraint numerically. As anticipated, 
the puncture tracks are sensitive to small changes in the 
initial data. Especially for three and more black holes 
changing the initial data, e.g. by solving the constraints 
rather than using an analytical approximation, can lead 
to qualitatively and quantitatively very different merger 
sequences. In any case, we confirmed the result of [30j, 
that, as expected, the puncture method lends itself 
naturally to the simulation of multiple black holes. 

Multiple black hole evolutions in numerical relativ- 
ity demand highly accurate methods for initial data and 
evolution schemes. Even in Newtonian dynamics, codes 
which calculate orbits require sophisticated methods to 
maintain stable and accurate evolution over long time 
scales. Evolutions of multiple black holes could be useful 
as a test case that taxes the accuracy of codes in com- 
paratively short evolutions. 

Simulations of three, four, or even more black holes 
lead to the following question about more general merger 
situations: How can we determine the number of black 
holes involved in a merger from the observation of grav- 
itational waves? A first analysis of this topic was given 
previously using a Newtonian approach [76,], with the 
perhaps surprising result that there are certain degen- 
eracies in the gravitational waves that prevent a trivial 
answer to this question. In the future, we plan to extend 
our research to a systematic study of the waveforms of 
multiple black hole configurations. 
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Appendix A: Convergence of high-order finite 
difference schemes for C" functions. 



In certain cases, the order of convergence of a finite 
difference scheme can be higher in the interior than at 
the boundary, without the lower order at the boundary 
spoiling the convergence in the interior (e.g. j77i |. Sec. 
2.12). Here we estimate the order of convergence of a 
standard p-order finite difference scheme for an elliptic 
problem, where the solution is C°° everywhere except on 
the origin where it is C" (where n < p). In order to 
simplify the notation, we will later restrict the examples 
to the one dimensional case. However, the extension to 
the three dimensional case is straightforward. 

Let C be an elliptic operator, i7 C an open domain, 
and u : 17 -> R the solution of the problem 



Cu{x) = p 
Bu{x) — Ui,{x) 



for 
for 



(Al) 
(A2) 



where S is a boundary operator, p : — > M is a source 
term, and u 6 C^{n) n C" (0). Let be a finite differ- 
ence representation of order p of £ in a mesh Vl^ C 
with a uniform grid size h. The numerical solution 
L/'' : 17'' ^ R satisfies 

C^U^{x^) ^ p^{x^) for f'^ell'', (A3) 
B^U^{x^)^u'i{x^) for x^cdn^, (A4) 

where is a discrete boundary operator and is the 
restriction of p on 17''. 

Given a point a; e 17, we identify points between 17'' 
and 17 by Xi — XQ + ih^ where i G {0, 1, . . . , N}. For every 
grid function we use as notation Ui := U'^{xi). The finite 
difference representation of C on the lattice location Xi 
has for each direction the form 



i+p 

E 

I—i—p 



ai~iUi, 



(A5) 



where the coefficients aj-i depend of the order of ap- 
proximation and the kind of stencil. For example, the 
standard 2nd order centered approximation to the sec- 
ond derivative is defined by oq = —2/h?, a±i = l/h?. 
The truncation error is defined by 



p 



(A6) 



where it'' is the restriction of u to the grid 17''. The 
approximation has the order of consistency p > if there 
is /iQ > which for all positive h < ho satisfies 



t'' < ChP, 



(A7) 



with a constant C > independent of h. The standard 
approach to analyzing the error in a finite difference ap- 
proximation is to expand each of the function values of 
m'' in a Taylor series about the point (xi). Taylor's the- 
orem states that for a function u G C"~^([a;i, a;]) and 
wGC'^((a:„x)), 



^-l (fc) 



u('^)(e) 



fe=0 



{x - X,)", (A8) 



where ^ G [xi,x] and u^"^^ denotes the n-th derivative. 
For grid functions the expansion formula is 



' ^ fc! 



n 



(A9) 



Using (jASP and (|A9p . it is possible to calculate 



i+p n— 1 
i+p n— 1 



,(fc) 



E E 

I—i—p k—0 



ai-i- 



(AlO) 



If n > p and the operator C contains a linear combination 
of derivatives up to order n — 1, then it is possible to 
select the coefficients at to cancel the remaining factors. 
We obtain 



i+p n— 1 



,(fe) 



I—i—p k—p 



i-\-p n~l 

E E«- 

I=i~p k=0 



u(")(C) 



(All) 



(j - «)"ft", 



where now the second summand starts at fc = p. If 
u(")(^)| ig bounded, the dominant term is of order h^. 
A substitution with (IA6I) leads to 



i+p (p) 

< I V aj^,\{I-i)P\hP, (A12) 
i=i-p ^ 

where the factor is bounded and independent of h. If we 
use the same scheme close to the origin, where n < p, we 
are not able to cancel terms lower than h": 



i+p n—1 



u(")(0 



C'v^=Ci4+ J2 E«^-'^tF(J-0"'^"- (A13) 



I=i-p k=0 
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The truncation error in this case is of order n < p, 



i+p 



I—i~p 

For example, for the operator 



(A14) 



(A15) 



the 4th order centered approximation to the second 
derivative is 



h„,h 



(A16) 



If w G C^(K)nC2(0) and e [a:i+i,a;i+2], a substit ution 
of the Taylor series of Xi±2 and Xi±i in equation (IA16|) 
results in 



1, /a^ h a3yh(^) 

h' 



dx^ 9 V 9a;3 



1 ^2 d'^H0 
18 



(A17) 



where we expand the term Xi+2 only up to 0{h'^). The 
truncation error is of order 0{h) . 
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